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I. INTRODUCTION 

Phenomena of synchronization of oscillators are wide 
spread in nature and technology. A variety of examples 
can be found in electronics, laser physics, biophysics, 
chemistry, neuroscience, etc. M- Synchronization in 
modern experiments for optomechanical, micromechani¬ 
cal, electronic oscillators is investigated, e.g. see [MU. 
The general picture of oscillatory modes in arrays of el¬ 
ementary oscillators substantially depends on a number 
of the elements and of the coupling type. In the simplest 
case of two oscillators with dissipative coupling one can 
observe mutual mode locking of the oscillators with dif¬ 
ferent natural frequencies, two-frequency quasi-periodic 
oscillations, the effect of ’’the oscillation death”, and the 
regime called the ’’broadband synchronization” specific 
for the oscillatory elements with non-identical control pa¬ 
rameters [ 1 H 3 , Il2l l6l|. With increase of a number of the 
oscillatory elements the picture becomes more compli¬ 
cated; many features of it have been established and un¬ 
derstood recently [rzH3. 

More complex is a case of reactive coupling (termed 
sometimes as the conservative coupling) pj-y, fl3hll5| . In 
radio-engineering and electronics, the coupling of such 
kind occurs in the case of presence of a reactive element 
(inductance) in the coupling circuit instead of a resistor 
giving rise to the dissipative coupling Q. A topical ex¬ 
ample of a system with the reactive coupling corresponds 
to ionic traps 0. In such traps ions are confined us¬ 
ing variable microwave fields, which restrict a magnitude 
of radial oscillations of the ions, and a constant electric 
field, limiting the axial motions. In a trap with many 
electrodes, the ions form a chain being located in the 
potential wells, and the nonlinearity provides the anhar- 
monic nature of oscillations of the ions in the wells. Ad¬ 
ditionally, the ions are irradiated by laser beams. The 
blue laser light of frequency larger than the natural fre¬ 
quency of the ion oscillations provides instability in the 


system. The red laser light of frequency less than the os¬ 
cillation frequency gives rise to dissipation. The ions in 
the chain are coupled due to the Coulomb repulsion. The 
simplest model of such a system is a chain of reactively 
coupled van der Pol oscillators [24|. In other concern the 
models of coupled van der Pol oscillators are studied in 
application to biological circadian rhythms [|25j and for 
arrays of nanoscale mechanical resonators |26j |. 

The reactive coupling is a phenomenon essentially 
more subtle than the dissipative coupling. The reason is 
that when constructing a reduced phase model one must 
take into account the second order effects in the coupling 
magnitude. In the case of two oscillators in the first or¬ 
der approximation the coupling effects are generally not 
manifested. For three or more oscillators the linear terms 
are present, but the resulting abridged system is conser¬ 
vative [271 l28j. 

The case of two oscillators was discussed in Refs. [IT 
l 15]. There the appropriate model for dynamics of the 
phase variable was derived and it was shown that the syn¬ 
chronization effect appears only with taking into account 
the terms of the second order in the coupling parameter. 
One more feature of the dynamics with the reactive cou¬ 
pling is the phase bistability; it means that depending 
on initial conditions the oscillators may synchronize ei¬ 
ther in phase, or in counter-phase. In the present Letter 
we consider dynamics and synchronization in a chain of 
three reactively coupled van der Pol oscillators. In con¬ 
trast to [23, [28|], we will assume non-identical oscillators 
for frequency and focus on the discussion of structure of 
parameter plane of frequency detuning. The first task is 
derivation of correct phase equations accounting all rel¬ 
evant effects. Then, they are studied using approaches 
developed in Refs. [22|, [23|. We reveal possible modes of 
complete and partial synchronization of the oscillators, 
the bifurcation mechanisms of destruction of the syn¬ 
chronization, and describe arrangement of the parameter 
domains of quasi-periodic modes with different number 
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of incommensurable frequencies. One of the main ques¬ 
tions we discuss concerns features and distinctions of the 
reactive coupling in comparison with the earlier known 
results for the dissipative coupling. 


II. THE PHASE MODEL 

Consider a set of equations 

x — (A — x 2 )x + x + e{x — y) = 0, 

V - (A - y 2 )y + (1 + Ai)y + e(y - x) + s(y - z) = 0, 

z — (A — z 2 )z + (1 + A 2 )z + e(z - y) = 0. 

(1) 

where A is a parameter controlling intensity of the self¬ 
oscillations, Ai and A 2 are the detuning parameters for 
the second and the third oscillators, and 5 is the coupling 
constant. 

If the excitation parameter A is small enough, as well 
as the detuning parameters, one can apply the slow- 
amplitude method for the analysis of the equations 0- 
Using the standard approach [1, 2]. one can derive the 
following equations for the real amplitudes ri and phases 
of the oscillators ^ (Landau-Stuart equations): 

2ri = 7*1 — r\ — er 2 sin#, 

2 r 2 = r 2 — 7*2 + et\ sin 9 — £r% sin p, 

2 r 3 = 7*3 — rf + er 2 sin p, 

2^1 = 5 — £^ COS#, (^) 

2 ip 2 = 26 + Ai - £^ cos # - cos p, 

2^3 =5 + A 2 -£^cosp. 

Here # = ^ 1 — ^ 2 , p = ^2 — 'ips are relative phases of the 
oscillators, and the parameters are normalized in respect 
to the small parameter A Hi H: 

r —^ \/Xr, t ^ tf A, 6 —)■ Ae, A —^ AA. (3) 

Subtracting pairwise the phase equations (2), we obtain 
the following equations for the relative phases: 

2# = —£ — Ai +e — ^) cosO + e^ cosip, 

2(f = £ + Ai - A 2 + e - ^) cos cp cos 6*. 

Let us set r* = 1 + f*, where the tilde designates per¬ 
turbations of the stationary orbits r = 1. The amplitude 
equations are strongly damped , so the orbits 

in a short time reach roughly the stationary amplitudes 
with some perturbations easily estimated from 0 as 

2fi = — 5 sin#, 2f2 = £ sin # — 5 sin p, 2r^ = £ sin ip. (5) 

In turn, from the phase equations 0 we obtain 

26 = —e — Ai + 2e(fi — V 2 ) cos # 

+ e(l + r 3 -f 2 )cos(^, 

2p = e + Ai - A 2 + 2 e{t 2 - rs) cos p 

— 5(1 + f \ — V 2 ) cos #. 


Substituting the expressions for the perturbations from 
(j5j) we get 

2# = — 5 — Ai + £ cos p — £ 2 sin 20 

+ £ 2 (sin p cos # — \ sin # cos p + \ sin 2 p) , 

2p = 6 + Ai — A 2 — ^ cos # — £ 2 sin 2 p 

+ 5 2 (sin # cos (/? — | sin p cos # + \ sin 2#) . 

(7) 

These are the correct phase equations for three reactively 
coupled oscillators derived up to the terms of order e 2 . 
Their structure is notably more complex than that for 
the dissipative coupling p], |22| . 

In contrast to the case of two oscillators [IM3 , the 
phase equations © do contain terms of the first order 
in the coupling strength £, however, for proper descrip¬ 
tion of the synchronization effects these terms are not 
sufficient. Indeed, if one neglect the quadratic terms, the 
matrix for perturbations of the stationary state of © is 


M = 


0, —esin p 

£ sin#, 0 


(8) 


The trace of this matrix is zero, S = 0. It means that 
a kind of “neutral” state occurs on the border between 
the stable and unstable solutions (conservative dynam¬ 
ics). Hence, in the system there is no main resonance at 
all in this approximation. As follows, for description of 
synchronization phenomena one has necessarily take into 
account the effects of the second order in the coupling 
parameter. 

Figure la illustrates bifurcations of equilibrium points 
of the system (0. In the case of dissipative coupling the 
border of the domain of complete synchronization corre¬ 
sponds to a curve of degenerate saddle-node bifurcations, 
where simultaneous merging occurs for a pair of the sad¬ 
dles with the stable and the unstable nodes [H, |29|. For 
the reactive coupling the curves of the saddle-node bifur¬ 
cations SN for merging of a stable node and a saddle and 
for an unstable node and a saddle do not coincide (33j| . 
Moreover, here a bifurcation of Andronov-Hopf is possi¬ 
ble (designated by H), where the equilibrium point be¬ 
comes unstable with appearance of the stable limit cycle 
departing from it. Therefore, the region of complete syn¬ 
chronization in the case of reactive coupling of phase os¬ 
cillators appears to be bounded both by the curves of the 
saddle-node bifurcations and of the line of the Andronov- 
Hopf bifurcations. Also we indicate in Figure la the 
points of codimension two: the cusp point CP and the 
Bogdanov-Takens point BT. In Figure lb, the domains 
inside of which system has one or two stable equilibria are 
shown using different colors. Thus, there is the simplest 
multistability. 

Figure 2 shows the chart of Lyapunov exponents [H, 
|23| of the system © on the parameter plane of frequency 
detuning of the oscillators (Ai,A 2 ). To draw the chart 
we compute two Lyapunov exponents of the system © 
Ai and A 2 at each pixel of the picture and attribute it 
with a color depending on the signature of the Lyapunov 
spectrum to visualize the following regimes: 
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FIG. 1: Bifurcation curves and points of the system 0, e = 
0.2. Digits in fragment b) indicate the number of coexisting 
stable equilibriums. 


a) Ai < 0, A 2 < 0 - the complete synchronization of 
three oscillators P(red), 

b) Ai = 0, A 2 < 0 - the two-frequency quasi-periodicity 
T 2 (yellow), 

c) Ai = 0, A 2 = 0 the three-frequency quasi-periodicity 
T 3 (blue). 

(Here the types of regimes are responsible to original sys- 
tem 0 . In this case it is convenient to compare the re¬ 
sults obtained at investigations of the phase model and 
the original system, see Fig.7.) 

The area of complete synchronization in Fig. 2 con¬ 
tains four “islands” [34j . In each island we observe a spe¬ 
cific kind of complete synchronization as illustrated in the 
phase portraits of Fig. 3 from (a) to (d). At the point (a) 
the relative phases are close to zero: 0 « 0, p « 0, and 
this is the synchronization mode of the in-phase type. 
At the point (b) the relative phases are 0 ~ —7r, p ~ 7 r; 
so, the first and the third oscillators are roughly in phase 
while the second oscillator is in the counter-phase rela¬ 
tively to them. This is the counter-phase synchroniza¬ 
tion. In the rest two islands we observe the complete 
synchronization of mixed type. In this case one of the 
pairs of the oscillators (1-2 or 2-3) are in phase, while 
the rest is in the counter-phase relative to them. The 
corresponding configurations of chain are shown in the 
right part of Fig.2|35|. 

Beside the region of complete synchronization, on the 
Lyapunov chart of Fig. 2 one can see a set of bands of two- 
frequency quasi-periodic regimes immersed in the domain 
of three-frequency regimes. Within each such band in¬ 
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FIG. 2: Chart of Lyapunov exponents on the parameter plane 
and configurations of basic modes of complete synchronization 
of the system 0 at e = 0.2. The letters designate points cor¬ 
responding to the phase portraits of Fig. 3. Used hereinafter 
color palette for the Lyapunov charts is shown below. 


variant curves of different types occur in the phase plane. 
Say, in the diagram 3(e) the relative phase of the first and 
the second oscillators 9 fluctuates around a certain equi¬ 
librium value, while the phase cj) varies across the whole 
range of values. This is a partial mode-locking of the first 
and the second oscillators. In the diagram 3(f) one ob¬ 
serves bounded oscillations of the relative phase of the 
second and thirds oscillators </>, so this is a partial mode¬ 
locking of the second and the third oscillators. 

We can classify regions of two-frequency modes with 
help of the rotation number w = p : q. Here p and q 
are numbers of intersection of the corresponding invari¬ 
ant curve with vertical and horizontal sides of the phase 
square. Only significant intersections should be used tak¬ 
ing into consideration 27r-periodicity of phase. So for 
Fig.3e the rotation number w = 0:1 (for both curves) 
and for Fig.3f w = 1:0. 

This classification becomes obvious if we compute a 
“torus map” using the numerical calculation of the fac¬ 
tors p and q at each point in the parameter plane. The 
corresponding illustration is given in Fig.4. We use the 
following rule of coloring. Blue color associates with 
regime of the rotation number w = 1:0. Decreasing 
of the rotation number corresponds to a piecemeal trans¬ 
formation of this color to green mode w = 1:1. Then 
green color is gradually transformed into the red for the 
mode with the rotation number w = 0:1. Stable equi¬ 
libriums are shown in white and the other modes - in 
black. Light gray color corresponds to the regime with 
contractible limit cycles when invariance curve has no 
significant intersections with the sides of the phase of a 
square. 
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FIG. 3: Phase portraits for the system of three reactively 
coupled oscillators © at £ = 0.2: a) Ai = 0, A 2 = 0; b) 
Ai = -0.45, A 2 = 0; c) Ai = -0.38, A 2 = -0.33; d) Ai = 0, 
A 2 = 0.36; e) Ai = -0.2, A 2 = 0.4; f) Ai = 0.05, A 2 = 0.22; 
g) Ai *= -0.2, A 2 = 0; h) Ai « -0.3, A 2 = 0. 



FIG. 4: Map of the two-frequency quasi-periodicity area of 
the system © 


Origin of the two frequency bands with a common sym¬ 


metry center in Figure 2 may be explained by the pres¬ 
ence of various possible resonances in the system. To 
substantiate this point, outline the partial frequencies uji 
of the system ([6|) . If we ” switch off’ in each equation all 
other oscillators, in the linear approximation we obtain 
from ([6]) 


£ , Ai s A 2 

wi =-, u) 2 =e+—, u) 3 = - + —- (9) 

A feature of the reactive coupling is that it shifts the 
partial frequencies by a value of order 5 (the shift for the 
central oscillator is twice as large as for the other ones 
because it interacts with two rather than one neighbors). 

Now, let us write out the main resonances for the par¬ 
tial frequencies. In round brackets the respective con¬ 
ditions are indicated in terms of the frequency detuning 
parameters: 


Wi = UJ 2 (Ai = -e ), 

w 2 = lo 3 (A 2 = Ai +e), , irV . 

<*>i = W3 (A 2 = 0), [W) 

u}\ + W 3 = 2tc>2 (A 2 = 2 Ai + 2 e) . 

The resonance conditions m determine the centers 
for the wide bands of two-frequency modes in Fig. 2 . Ad¬ 
ditional resonances are possible too, say, u)\ + UJ 2 = 2 cj 3 
etc. They correspond to more narrow bands in Fig. 2. 
The resonances of higher order give rise to invariant 
curves with larger number of intersections with the sides 
of the phase square. 

Note that the resonances uj\ = uj 3 and uj\ T CJ 3 = 2 cj 2 
determine lines of symmetry on the parameter plane, see 
Figs. 1 and 2 . This is due to symmetry of these resonance 
conditions in respect to permutation of the first and third 
oscillators. Say, under the condition uj 1 = CJ 3 , i.e. A 2 = 0 , 
the equations © are transformed one to other under the 
variable change 9 —ip. As a result, the phase portraits 
are symmetric about the line p = —0 , see Fig. 3g,h. 
Analogously, with the condition A 2 = 2 (Ai + e) the 
equations are invariant in respect to the variable change 

6 p + 7r. 

A case is interesting of equality of all the partial fre¬ 
quencies (jJ\ — UJ 2 = ^3 that corresponds to the sym¬ 
metry center on the parameter plane Ai = —£^2 = 0. 
The phase portrait at this point is shown in Fig. 3g, 
and at the point close to it at - in Fig. 3h. In this case 
invariant curves of other kind arise, which are distinct 
in their topological properties. The curves in Fig. 3g 
may be called contractible , while those in diagrams (a)- 
(/) are called rotational [3l| (for the dissipative cou¬ 
pling the second case is typical [22|). In the first case 
we have a limit cycle going around the unstable equi¬ 
librium point. Then, both the relative phases fluctuate 
about some mean value. This mode may be characterized 
as a partial mode-locking of all three oscillators. The fre¬ 
quency spectrum produced by the system (pQ) will contain 
not only the basic frequency, but also a set of components 
associated with an additional new time scale, the period 
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of travel of the representative point around the limit cycle 
of the phase model. 

Contractible limit cycles can occur, as we have already 
noted, as a result of Andronov-Hopf bifurcation. In ad¬ 
dition, the system can demonstrate nonlocal bifurcations 
when the resulting limit cycle arises from the separatrix 
loop of the saddle. 

Note that several types of multistability are possi¬ 
ble in the system with reactive coupling as seen from 
Fig. 3. Diagrams (a)-(d) correspond to situation when 
an equilibrium state in the phase space corresponding to 
the complete synchronization coexists with an invariant 
curve corresponding to a two-frequency quasi-periodic 
regime [36]. In diagram (e) two regimes coexist (the in- 
phase mode and the counter-phase mode) corresponding 
to partial mode-locking of the first and second oscilla¬ 
tors. In turn, in diagram (g) there coexist two regimes 
of partial synchronization of all three oscillators. 

Multistability affects the appearance of the Lyapunov 
chart. Namely, when choosing different initial conditions 
one can observe either periodic or quasi-periodic regimes, 
as can be seen from a comparison of Figures 1 and 2. 


III. DYNAMICS OF THE ORIGINAL SYSTEM 


0.05 




FIG. 5: Bifurcation curves and points of the system (1), A = 
0.1, £ = 0.02. 


Let us illustrate the effectiveness of the phase model. 
For this purpose we represent the results of the bifurca¬ 
tion analysis of the original system © at A = 0.1, Figure 
5. In accordance with the renormalization rules (3) the 
coupling parameter £ = 0.02 is selected, which allows to 
compare Figure 5 for the original system and Figure 1 
for the phase model. Now instead of a saddle-node bifur¬ 
cation of equilibria there is a corresponding bifurcation 
of limit cycles SNC, instead of the Andronov-Hopf bi¬ 
furcation - Neimark-Sacker bifurcation NS, and instead 
of Bogdanov-Takens points - points of 1:1 resonance Rl. 
From Figure 1 and Figure 5 we can see that the Figures 
are similar to each other, which indicates the effectiveness 
of the phase model for these values of coupling parame¬ 
ter. 

It should be noted that the effectiveness of the phase 
model will increase with a decrease in the coupling pa¬ 
rameter e. With increasing of coupling parameter its effi¬ 
ciency falls. Fig.6 illustrates this fact, demonstrating bi¬ 
furcation lines of the phase model and the original system 
©• For the phase model we select coupling parameter 
£ = 0.6, so for original system we have £ = 0.06 (taking 
into account A = 0.1 rescaled on ©)• Some common 
features - the presence of four lobes - remain. However, 
the pictures are different in details. Namely, external 
boundaries of lobes now are mainly the Neimark-Saker 
lines NS. Lines of saddle-node bifurcations of limit cy¬ 
cles SNC form only small segments of the boundary of 
complete synchronization area in the vicinity of the cusp 
points CP associated with bistability areas. One more 
new feature is the appearance of a double Neimark-Sacker 
bifurcation points NS-NS. 



-0.26 Ai 0.14 


FIG. 6: a) Bifurcation lines and dots of phase model ©, 
e = 0.6; b) the similar illustration for the original system ([I]). 
A = 0.1, £ = 0.06. 


Chart of Lyapunov exponents of the original system © 
is presented in Fig.7 for A = 0.1, £ = 0.06. An enlarged 
fragment of this chart in Fig.7b should be compared with 
Fig.6b. 
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-0.6 . A, 0.6 

a) 


- 0.6 



FIG. 7: The chart of Lyapunov exponents and its magnified 
fragment of the system (pQ), A = 0.1, £ — 0.06. 

With the growth of the A control parameter both the 
Landau-Stuart and phase model will work bad to worse. 
In fig.8a the Lyapunov chart is shown for the system 
© for the case A=l, 5=0.6. Now the picture is much 
more complex. Fig.8b shows a new effect: at the inter¬ 
section of numerous bands of two-frequency modes there 
are located regions of higher resonances, to which there 
correspond periodic regimes. This is a characteristic res¬ 
onance Arnold web J32[. Letter D denotes the area of the 
trajectories’ escape to infinity. 

IV. CONCLUSION 

The results presented here may be of interest for sys¬ 
tems such as ion traps biophysical systems (25) . 

etc. At the same time, due to universality of our mod¬ 
els, they are important as well in the general theory of 
synchronization. The main result is that in the descrip¬ 
tion of reactive coupling in the framework of the slow 
amplitudes approach it is of principal importance to ac¬ 
count the effects of the second order in the coupling con¬ 
stant. The region of complete synchronization consists 
of four ” islands” in which the synchronization modes are 
observed corresponding to in-phase, counter-phase, and 


mixed types of oscillations in the chain. The bifurcation 



FIG. 8: The chart of Lyapunov exponents and its magnified 
fragment of the system ©, A=l, 5 = 0.6. 


picture of the model formulated in terms of phase vari¬ 
ables is significantly different from the case of dissipative 
coupling. In particular, Andronov - Hopf bifurcation is 
possible responsible for occurrence of partial synchroniza¬ 
tion regimes of all three oscillators. For the arrangement 
of the parameter plane of the frequency detunings, res¬ 
onances between the partial frequencies are important, 
and for the reactive coupling the resonance conditions 
appear to depend on the coupling strength. One more 
feature is that the quasi-periodic modes of various types 
can co-exist with complete synchronization, giving rise 
to many kinds of multistability (coexistence of different 
attractors). 
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